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Optical imaging relies on the ability to illuminate an object, collect and analyze the light it scatters 
or transmits. Propagation through complex media such as biological tissues was so far believed to 
degrade the attainable depth as well as the resolution for imaging because of multiple scattering. This 
is why such media are usually considered opaque. Very recently, we have proven that it is possible to 
measure the complex mesoscopic optical transmission channels that allows light to traverse through 
such an opaque medium. Here we show that we can optimally exploit those channels to coherently 
Qh transmit and recover with a high fidelity an arbitrary image, independently of the complexity of the 
propagation. 

In a classical optical system, the propagation of a complex field from one plane to another is well understood, be it 
by Fresnel or Fraunhofer diffraction theory, or ray-tracing for more complex cases [1 j. However, all these approaches 
break down when multiple scattering occurs [2 . A medium where light is scattered many times mixes in a seemingly 
random way all input k- vectors, and is usually considered opaque. Until recently, scattering has always been considered 
as noise |3j and most imaging techniques in turbid media rely on ballistic photons only [4, 5 which prevents the study 
of thick scattering samples. Following works in acoustics [6], recent experiments have demonstrated that multiply 
scattered light can nonetheless be harnessed, thanks to wavefront control [7H2], and even put to profit to surpass what 
one can achieve within a homogenous medium in terms of focusing [10] . 

In our experiment (see figure HI), we illuminate with a laser an object (displayed via a spatial light modulator, 
or SLM), and recover its image on a CCD camera, after propagation through a thick opaque sample. As expected, 
we measure on the camera a speckle, that bears no resemblance to the original image. This speckle is the result of 
multiple scattering and interferences in the sample. 

Although it can be described on average by diffusion equation or Monte Carlo simulation [11], the propagation 
O ! through a real linear multiple scattering medium is too complex to be described by classical means. Nonetheless, 
multiple scattering is deterministic and information is not lost. In other terms, the measured pattern on the CCD is 
the result of the transmission of light through a large number of very complicated optical channels, each of them with 
a given complex transmission. Here, we study the inverse problem of the reconstruction of an arbitrary image, and 
show that it is possible to recover it through the opaque medium. A prerequisite is however to measure the so-called 
CO transmission matrix (TM) of our optical system. 

We define the mesoscopic transmission matrix (TM) of an optical system for a given wavelength as the matrix K of 
the complex coefficients k mn connecting the optical field (in amplitude and phase) in the m th of M output free mode 
to the one in the n th of N input free mode. Thus, the projection E™ 1 of the outgoing optical field on the m th free 
mode is given by E™ 1 = J2 n kmnE™ where E™ is the complex amplitude of the optical field in the n th incoming free 
mode. In essence, the TM gives the relationship between input and output pixels, notwithstanding the complexity of 
the propagation, as long as the medium is stable. A Singular Value Decomposition (SVD) of the TM gives the input 
and output eigenmodes of the system and singular values are the amplitude transmission of these modes. 

Inspired by various works in acoustics [I2j [13] and electromagnetism [14] , we demonstrated in [15] that it is possible 
to measure the TM of a linear optical system that comprises a multiple scattering medium. In a nutshell, we send 
several different wavefront s with the SLM, record the results on the CCD, and deduce the TM using phase-shifting 
interferometry. The singular value distribution of a TM of a homogeneous zone of the sample follows the quarter-circle 
law (i.e. there is no peculiar input/output correlation [16]) which indicates that light propagation is in the multiple 
scattering regime with virtually no ballistic photons left. 

Using this technique, we have access to if 6s = K x S re f, where S re f is a diagonal matrix due to a static reference 
speckle. The input and output modes are the SLM and the CCD pixels respectively. The measured matrix Ko^ is 
sufficient to recover an input image. This TM measurement takes a few minutes, and the system is stationary well 
over this time. Once the matrix is measured, we generate an amplitude object Eq^j by subtracting two phase objects 
(see Methods for details). A realization takes a few hundred ms, limited only by the speed of the SLM. 

Here, our aim is to use the TM to reconstruct an arbitrary image through the scattering sample: we need to 
estimate the initial input E ^ from the output amplitude speckle E out . This problem consists in using an appropriate 
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Fig. 1. Experimental setup. A 532 nm laser is expanded and reflected off a spatial light modulator (SLM) . 
The laser beam is phase-modulated then focused on the multiple-scattering sample and the output intensity speckle 
pattern is imaged by a CCD-camera. L, lens. P, polarizer. D, diaphragm. The object to image is synthetized directly 
by the SLM, and reconstructed from the output speckle thanks to the transmission matrix. 

combination of the medium channels and therefore using a weighting of singular modes/singular values of the TM 
matched to the noise and to the transmitted image. Noises of different origins (laser fluctuation, CCD readout noise, 
residual amplitude modulation) degrade the fidelity of the TM measurement. It is the exact analog of Multiple-Input 
Multiple-Output (MIMO) information transmission in complex environment that have been studied in the past few 
years in wireless communications [17]. This inverse problem also bears some similarities to optical tomography [HJ[T9], 
albeit in a coherent regime [20] . 

There are two straightforward options, (i) Without noise, a perfect image transmission can be performed by the 
use of the inverse matrix (or pseudo-inverse matrix for any input /output pixels ratio) since K~^ s K b s = / where I is 
the identity matrix. Unfortunately, this operator is very unstable in presence of noise. Singular values of K~^ s are the 
inverse ones of if 6s> thus singular values of K & s below noise level result in strong and aberrant contributions. The 
reconstructed image can hence be unrelated with the input, (ii) In a general case, another possible operator for image 
transmission is the Time Reversal operator. This operator is known to be stable regarding noise level since it takes 
advantage of the strong singular values to maximize energy transmission [12] . Its monochromatic counterpart is phase 
conjugation (classically used to compensate dispersion in optics [21]) and is performed using K* ohs . Kl bs K b s has a 
strong diagonal but the rest of it is not null, which implies that the fidelity of the reconstruction rapidly decreases 
with the complexity of the image to transmit [22 . A more general approach is to use a Mean Square Optimized 
operator (MSO), that we note W. This operator minimizes transmission errors [13 [24], estimated by the expected 

value E \ [W.E out — E tn ~\ [W.E out — E tn ~\ >. For an experimental noise of standard deviation a on the output pixels, 
W reads : 



W = 



K bs -Ko bs 



a.I 



K 



obs' 



(i) 



Without noise, W reduces to the inverse matrix K~^~, which is optimal in this configuration, while for a very high 
noise level it becomes proportional to the transpose conjugate matrix K* ohs , the phase conjugation operator. It is 
important to note that a has the same dimension as K] )hs .K } )S and thus has to be compared with the square of 
the singular values of K b s . Because of this experimental noise the reconstruction is imperfect. We estimate the 
reconstruction fidelity by computing the correlation between the image and the object. 

A general principle is that the reconstruction noise can be lowered by increasing the number of degree of freedom 
(Njjqf) that we measure and control. For a given object corresponding to N input pixels, we investigated two 
possibilities: averaging over disorder realizations and increasing the number of output modes M. 

A possible way to average over disorder is to illuminate the object with different wavefronts. It is formally equivalent 



to transmitting the same image through different channels as if the image propagated through different realizations of 
disorder. To that end, we use different combinations of random phase masks to generate the same 'virtual object' (see 
methods). We use this technique to virtually increase N^of-, and we average the results to lower the reconstruction 
noise. It is the mochromatic equivalent of using broadband signals, which takes advantage of temporal degree of 
freedom [25]. We show in figure 2 the results for the image transmission of a gray-scale 32 by 32 pixels pattern, 
and detected on a 32 by 32 pixels region on the CCD. We tested MSO at different noise level for one realization 
and for averaging over 40 'virtual realizations' with random phase masks. To find the optimal MSO operator, we 
numerically compute the optimal a that maximizes the image reconstruction, hence obtaining an estimation of the 
experimental noise level. A simple inverse filtering does not allow image reconstruction, even with averaging, while 
phase conjugation converges toward 75% correlated image. In contrast, optimal MSO, allowed a 94% correlation for 
40 averaging (and a modest 34% correlation in one realization). In addition, Optimal MSO is very robust to the 
presence of ballistic contributions that strongly hinder reconstruction in phase conjugation (see discussion). 
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Fig. 2. Comparisons of the reconstruction methods, a. initial gray scale object and b. a typical output speckle 
figure after the opaque medium, c. and f. are experimental images obtained with one realization using respectively 
phase conjugation and MSO operator, d., e. and g. are experimental images averaging over 40 'virtual realizations' 
using respectively inverse matrix, phase conjugation and MSO operator. Values in insets are the correlation with the 
object a., h. Correlation coefficient between E irng and E Q ^ as a function of ^fo (line) and singular value distribution 
of Kobs (bars). Results are obtained averaging over 100 'virtual realizations' of disorder and both ^[a and singular 
values share the same scale on the abscissa axis. 



The second approach in order to add degrees of freedom is to increase the number M of independent pixels recorded 
on the CCD. In contrast with focusing experiments where the quality of the output image depends on the number of 
input modes N |5], the quality of image reconstruction depends on the number of output modes M. An important 
advantage is that the limiting time in our experiment is the number or steps required to measure the TM, equal to 4N. 
Thus, we can easily increase M by increasing the size of the image recorded without increasing the measurement time. 
More than just modifying the Ndof-, the ratio 7 = M/N > 1 is expected to change the statistics of the TM. Random 
Matrix Theory (RMT) predicts that for those matrices the smallest normalized singular value reads A^ = (1 — \A/t) 
[16j [23] . If we increase 7 so does the minimum singular value A^. In a simple physical picture, recording more 
information at the output results in picking between all available channels those that convey more energy through 
the medium. If the energy transported by the most inefficient channel reaches and exceeds the noise level, the TM 
recording is barely sensitive to the experimental noise. We expect that for an appropriate ratio 7, A^ reaches the 
experimental noise level. At this point, no singular values can be drowned in the noise and the pseudo-inverse operator 
can be efficiently used. 

We experimentally recorded the TM for different values of 7 > 1 and tested optimal MSO and pseudo-inversion. 
The results are shown in Figure [3j Adding degrees of freedom strongly improves the quality of the reconstruction. 
We see that the quality of the reconstructed image increases with 7 and reaches a > 85% fidelity for the largest 
value of 7 = 11, without any averaging. The minimum singular value A^ also increases with 7. As expected, for 
A^ > <j opt , pseudo-inversion is equivalent to optimal MSO. One notice that experimental A^ are always smaller than 
their theoretical predictions. This deviation can be explained by the amplitude of the reference pattern \S re f\ that 
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Fig. 3. Influence of the number of output detection modes, a. Correlation coefficient between Ei mg and 
Eobj as a function of the asymmetric ratio 7 = M/N of output to input pixels for MSO (dashed line) and for pseudo- 
inversion (solid line), without any averaging. Error bars correspond to the dispersion of the results over 10 realizations. 
b. Experimental (solid line) and Marcenko-Pastur [16] predictions (dashed line) for the minimum normalized singular 
value as function of 7. The horizontal line show the experimental noise level Vcr opt • 



induces correlations in the matrix. It is well known in RMT that correlations modify the SVD of a matrix of identically 
distributed elements [TS] . 

So far, we tested image transmission in the case of a homogeneous medium, but what would be the results in more 
complex conditions ? Here, we study the robustness of this technique in presence of ballistic contributions, that is, 
a fraction of light that have not been scattered at all. The singular values of if 6s are proportional to the amplitude 
transmitted through each channel of the system. Ballistic contributions should give rise to strong singular values 
corresponding to the apparition of channels of high transmission. These are not spatially homogeneously distributed 
in energy, contrarily to multiply scattered contributions. Phase conjugation maximizes energy transmission in channel 
of maximum transmission [12]. Therefore, ballistic high singular values contributions should be predominant in phase 
conjugation, independently the image Eq^ and will not efficiently contribute to image reconstruction. MSO should 
not be affected since it reaches the optimum intermediate between inversion which is stable except for singular values 
below noise level, and phase conjugation which forces energy in maximum singular value channels. In other words, 
MSO will lower weight of channels which do not efficiently contribute to the image reconstruction. 

To experimentally study this effect, we moved the collection objective closer to the sample on a thinner and 
less homogeneous region, where some ballistic light could be recorded. We study in Figure 4 the quality of the 
reconstruction as function of a for both experimental conditions (with and without ballistic contributions). Both 
experiments give comparable results with 93.6% and 94.5% correlation coefficient with the optimal MSO operator 
and both give very low correlation results for inverse matrix operator. With the phase conjugation operator (equivalent 
to MSO for high a), the experiment sensitive to ballistic contributions, give a low correlation coefficient around 35%, 
to be compared to the value of 75.7% we obtained through the multiple scattering sample. This difference can be 
explained by the presence of a few high singular values contributions (two times greater than the maximum of the 
other singular values) that perturbate the image reconstruction. 

To conclude, we have shown that the transmission matrix allows a rapid and accurate reconstruction of an arbitrary 
image after propagation through a strongly scattering medium (see supplementary movie). Our approach gives a 
general framework for coherent imaging in complex media, going well beyond focusing and phase conjugation. It is 
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Fig. 4. Influence of the transmission channels on the reconstruction Correlation coefficient between Ei mg and 
Eobj as a function of y/a (line) and singular value distribution of K i, s with ballistic contributions in the transmission 
matrix, averaged over 100 'virtual realizations' of disorder. Both y/a and singular values share the same scale on 
the abscissa axis. Ballistic contributions strongly degrade the reconstruction in phase conjugation while MSO is 
unaffected. 

valid for any linear complex media, and could be extended to several novel photonic materials, whatever the amount of 
scattering or disorder (from complete disorder to weakly disordered photonic crystals [27] , and from super diffusive [28] 
to Anderson localization [29]). The quality of the reconstruction can be increased by harnessing the degrees of freedom 
of our system, and is very resilient to noise. In addition to its obvious interest for imaging, this experiment strikingly 
shows that manipulation of wave in complex media is far from limited to single or multi-point focusing. In particular, 
due to spatial reciprocity, a similar experiment could be performed using an amplitude and phase modulator by 
shaping the input wavefront to form an image at the output of an opaque medium, which would allow a resolution 
solely limited by the numerical aperture of the scattering medium [10] . The main current limitation is the speed of 
the TM measurement, which is limited only by the spatial light modulator. Nevertheless faster technologies emerge, 
such as micromirror arrays or ferromagnetic SLMs, that might in the future widen the range of application domains 
for this approach, including the field of biological imaging. 

Methods 



Imaging Setup 

The experimental setup consists on an incident light from a 532 nm laser source (Laser Quantum Torus) that is 
expanded, spatially modulated by a Spatial Light Modulator (Holoeye LC-R 2500), and focused on an opaque strongly 
scattering medium : 80 ± 25 jam thick deposit of ZnO ( Sigma- Aldrich 96479) with a measured transport mean free 
path of 6±2 microns on a standard microscope glass slide. Polarization optics select an almost phase-only modulation 
mode [26] of the incident beam, with less than 10% residual amplitude modulation. The surface of the SLM is imaged 
on the pupil of a lOx objective, thus a pixel of the SLM matches a wave vector at the entrance of the scattering 
medium. The beam is focused at one side of the sample and the output intensity speckle is imaged on the far side 
(0.3 mm from the surface of the sample) by a 40x objective onto a 10-bit CCD camera (AVT Dolphin F-145B). 

Generation of the amplitude object 

Since there is no simple way to control the amplitude and phase of the incident beam, we generate a virtual amplitude 
object (E ob i with sffl G [0,1]) by substracting two phase objects. This method is more flexible than placing a 
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real amplitude object in the plane of the SLM. From any phase mask E^ ase we could generate a second mask E y phase 
where the phase of the m th pixel is shifted by sffl .7r. We have e 
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